{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Untitled17.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.8"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/jeshraghian/snntorch/blob/master/examples/tutorial_2_neuronal_dynamics.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HzIQBw28NL8h"
      },
      "source": [
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/snntorch_alpha_w.png?raw=true' width=\"400\">\n",
        "\n",
        "# snnTorch - Neuronal Dynamics with ``snntorch``\n",
        "## Tutorial 2\n",
        "### By Jason K. Eshraghian (www.jasoneshraghian.com)\n",
        "\n",
        "<a href=\"https://colab.research.google.com/github/jeshraghian/snntorch/blob/tutorials/examples/tutorial_2_neurons.ipynb\">\n",
        "  <img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/>\n",
        "</a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ep_Qv7kzNOz6"
      },
      "source": [
        "# Introduction\n",
        "In this tutorial, you will:\n",
        "* Learn the fundamentals of the leaky integrate-and-fire (LIF) neuron model\n",
        "* Use snnTorch to implement variations of the LIF model: \n",
        "  * Lapicque's neuron model (1st order)\n",
        "  * Synaptic conductance-based model (2nd order)\n",
        "  * Alpha model (a hacked version of the Spike Response model)\n",
        "<!-- * Plot the output behavior of the neurons -->\n",
        "<!-- * Interpret the computational graph of a spiking neuron -->\n",
        "<!-- * Automatically initialize the hidden states of the neurons [keep in tute, but delete explanation]? -->\n",
        "* Implement a feedforward spiking neural network\n",
        "\n",
        ">Part of this tutorial was inspired by the book [*Neuronal Dynamics:\n",
        "From single neurons to networks and models of cognition*](https://neuronaldynamics.epfl.ch/index.html) by\n",
        "Wulfram Gerstner, Werner M. Kistler, Richard Naud and Liam Paninski.\n",
        "\n",
        "If running in Google Colab:\n",
        "* You may connect to GPU by checking `Runtime` > `Change runtime type` > `Hardware accelerator: GPU`\n",
        "* Next, install the latest PyPi distribution of snnTorch by clicking into the following cell and pressing `Shift+Enter`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SPQITvDuNNJg"
      },
      "source": [
        "!pip install snntorch"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xmTt5dvyXNNy"
      },
      "source": [
        "# 1. The Spectrum of Neuron Models\n",
        "A large variety of neuron models are out there, ranging from biophysically accurate models (i.e., the Hodgkin-Huxley models) to the extremely simple artificial neuron that pervades all facets of modern deep learning.\n",
        "\n",
        "**Hodgkin-Huxley Neuron Models**$-$While biophysical models can reproduce electrophysiological results with a high degree of accuracy, their complexity makes them difficult to use. We expect this to change as more rigorous theories of how neurons contribute to higher-order behaviors in the brain are uncovered.\n",
        "\n",
        "**Artificial Neuron Model**$-$On the other end of the spectrum is the artificial neuron. The inputs are multiplied by their corresponding weights and passed through an activation function. This simplification has enabled deep learning researchers to perform incredible feats in computer vision, natural language processing, and many other machine learning-domain tasks.\n",
        "\n",
        "**Leaky Integrate-and-Fire Neuron Models**$-$Somewhere in the middle of the divide lies the leaky integrate-and-fire (LIF) neuron model. It takes the sum of weighted inputs, much like the artificial neuron. But rather than passing it directly to an activation function, it will integrate the input over time with a leakage, much like an RC circuit. If the integrated value exceeds a threshold, then the LIF neuron will emit a voltage spike. The LIF neuron abstracts away the shape and profile of the output spike; it is simply treated as a discrete event. As a result, information is not stored within the spike, but rather the timing (or frequency) of spikes. Simple spiking neuron models have produced much insight into the neural code, memory, network dynamics, and more recently, deep learning. The LIF neuron sits in the sweet spot between biological plausibility and practicality. \n",
        "\n",
        "<center>\n",
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/examples/tutorial2/2_1_neuronmodels.png?raw=true' width=\"1000\">\n",
        "</center>\n",
        "\n",
        "<!-- Researchers might spend their entire lives dedicated to developing neuron models. Some of these models are straightforward extensions of the HH and LIF models, while other models find completely different applications, such as in neuropharmocology. -->\n",
        "\n",
        "The different versions of the LIF model each have their own dynamics and use-cases. snnTorch currently supports four types of LIF neurons:\n",
        "* Lapicque's RC model: ``snntorch.Lapicque``\n",
        "* Non-physical 1st order model: ``snntorch.Leaky`` \n",
        "* Synaptic Conductance-based neuron model: ``snntorch.Synaptic``\n",
        "* Alpha neuron Model: ``snntorch.Alpha``\n",
        "\n",
        "Before learning how to use them, let's understand how to construct a simple LIF neuron model.\n",
        "\n",
        "<!-- In general, the most obvious difference is that the SRM0 model incorporates a delay between the input and output. When an input spike arrives at an SRM0 neuron, the membrane potential will increase over a finite time. If an output spike were to be triggered, it would experience a delay with respect to the input. On the other hand, Stein's model allows for an instantaneous rise of membrane potential. We'll dig into where these might be useful shortly. -->\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Nea8oBorr_KZ"
      },
      "source": [
        "# 2. The Leaky Integrate-and-Fire Neuron Model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YKsrN5feQ2Dz"
      },
      "source": [
        "## 2.1 Spiking Neurons: Intuition"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YGStnfjzsGKb"
      },
      "source": [
        "A neuron might be connected to 1,000 $-$ 10,000 other neurons. If one neuron spikes, all of these downhill neurons will feel it. But what determines whether a neuron spikes in the first place? The past century of experiments demonstrate that if a neuron experiences *sufficient* stimulus at its input, then we might expect it to become excited and fire its own spike. \n",
        "\n",
        "Where does this stimulus come from? It could be from:\n",
        "* the sensory periphery, \n",
        "* an invasive electrode artificially stimulating the neuron, or in most cases,\n",
        "* from other pre-synaptic neurons. \n",
        "\n",
        "<center>\n",
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/examples/tutorial2/2_2_intuition.png?raw=true' width=\"500\">\n",
        "</center>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TaVTRojJ2Dl_"
      },
      "source": [
        "Given that these spikes are very short bursts of electrical activity, it is quite unlikely for all input spikes to arrive at the neuron body in precise unison. This indicates the presence of temporal dynamics that 'sustain' the input spikes, kind of like a delay.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "l44jI7A2ReB_"
      },
      "source": [
        "## 2.2 The Passive Membrane"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Od5tfXv6_wcq"
      },
      "source": [
        "Like all cells, a neuron is surrounded by a thin membrane. This membrane is a lipid bilayer that insulates the conductive saline solution within the neuron from the extracellular medium. Electrically, the two conductors separated by an insulator act as a capacitor. \n",
        "\n",
        "Another function of this membrane is to control what goes in and out of this cell (e.g., ions such as Na$^+$). The membrane is usually impermeable to ions which blocks them from entering and exiting the neuron body. But there are specific channels in the membrane that are triggered to open by injecting current into the neuron. This charge movement is electrically modelled by a resistor.\n",
        "\n",
        "<center>\n",
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/examples/tutorial2/2_3_passivemembrane.png?raw=true' width=\"450\">\n",
        "</center>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wT_S4b4HL2ir"
      },
      "source": [
        "Now say some arbitrary time-varying current $I_{\\rm in}(t)$ is injected into the neuron, be it via electrical stimulation or from other neurons. The total current in the circuit is conserved, so:\n",
        "\n",
        "$$I_{\\rm in}(t) = I_{R} + I_{C}$$\n",
        "\n",
        "From Ohm's Law, the membrane potential measured between the inside and outside of the neuron $U_{\\rm mem}$ is proportional to the current through the resistor:\n",
        "\n",
        "$$I_{R}(t) = \\frac{U_{\\rm mem}(t)}{R}$$\n",
        "\n",
        "The capacitance is a proportionality constant between the charge stored on the capacitor $Q$ and $U_{\\rm mem}(t)$:\n",
        "\n",
        "\n",
        "$$Q = CU_{\\rm mem}(t)$$\n",
        "\n",
        "The rate of change of charge gives the capacitive current:\n",
        "\n",
        "$$\\frac{dQ}{dt}=I_C(t) = C\\frac{dU_{\\rm mem}(t)}{dt}$$\n",
        "\n",
        "Therefore:\n",
        "\n",
        "$$I_{\\rm in}(t) = \\frac{U_{\\rm mem}(t)}{R} + C\\frac{dU_{\\rm mem}(t)}{dt}$$\n",
        "\n",
        "$$\\implies RC \\frac{dU_{\\rm mem}(t)}{dt} = -U_{\\rm mem}(t) + RI_{\\rm in}(t)$$\n",
        "\n",
        "The right hand side of the equation is of units **\\[Voltage]**. On the left hand side of the equation, the term $\\frac{dU_{\\rm mem}(t)}{dt}$ is of units **\\[Voltage/Time]**. To equate it to the left hand side (i.e., voltage), $RC$ must be of unit **\\[Time]**. We refer to $\\tau = RC$ as the time constant of the circuit:\n",
        "\n",
        "$$ \\tau \\frac{dU_{\\rm mem}(t)}{dt} = -U_{\\rm mem}(t) + RI_{\\rm in}(t)$$\n",
        "\n",
        "The passive membrane is therefore described by a linear differential equation.\n",
        "\n",
        "For a derivative of a function to be of the same form as the original function, i.e., $\\frac{dU_{\\rm mem}(t)}{dt} \\propto U_{\\rm mem}(t)$, this implies the solution is exponential with a time constant $\\tau$.\n",
        "\n",
        "Say the neuron starts at some value $U_{0}$ with no further input, i.e., $I_{\\rm in}(t)=0$. The solution of the linear differential equation is:\n",
        "\n",
        "$$U_{\\rm mem}(t) = U_0e^{-\\frac{t}{\\tau}}$$\n",
        "\n",
        "The general solution is shown below.\n",
        "\n",
        "<center>\n",
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/examples/tutorial2/2_4_RCmembrane.png?raw=true' width=\"500\">\n",
        "</center>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rTXS_aSmRs-3"
      },
      "source": [
        "## 2.3 Lapicque's LIF Neuron Model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Jt2q6tZiWkgT"
      },
      "source": [
        "This similarity between nerve membranes and RC circuits was observed by [Louis Lapicque in 1907](https://core.ac.uk/download/pdf/21172797.pdf). He stimulated the nerve fiber of a frog with a brief electrical pulse, and found that membranes could be approximated as a capacitor with a leakage. We pay homage to his findings by naming the basic LIF neuron model in snnTorch after him. \n",
        "\n",
        "Most of the concepts in Lapicque's model carry forward to other LIF neuron models. Now let's simulate this neuron using snnTorch."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WtveAGG0zE0n"
      },
      "source": [
        "### 2.3.1 Lapicque: Without Stimulus"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "um1s01gTzUC0"
      },
      "source": [
        "First, import the packages needed to run Lapicque's neuron model: snnTorch and PyTorch."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sUyODsBtWkAG"
      },
      "source": [
        "import snntorch as snn\n",
        "import torch"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6zpgNEbizd8s"
      },
      "source": [
        "The membrane potential has a time constant $\\tau = RC$ associated with it. This can be equivalently represented by a decay rate $\\beta$ that specifies the ratio of potential between subsequent time steps:\n",
        "\n",
        "$$\\beta = \\frac{U_0e^{-\\frac{1}{\\tau}}}{U_0e^{-\\frac{0}{\\tau}}} = \\frac{U_0e^{-\\frac{2}{\\tau}}}{U_0e^{-\\frac{1}{\\tau}}} = \\frac{U_0e^{-\\frac{3}{\\tau}}}{U_0e^{-\\frac{2}{\\tau}}}=~~...$$\n",
        "$$\\implies \\beta = e^{-\\frac{1}{\\tau}}$$\n",
        "\n",
        "Setting $\\tau = 5\\times 10^{-3} \\implies \\beta \\approx 0.819$:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ohshwOCU6Vbm"
      },
      "source": [
        "# RC time constant\n",
        "tau_mem = 5e-3\n",
        "time_step = 1e-3 # one time step = 1ms\n",
        "\n",
        "# decay p/time step\n",
        "beta = float(torch.exp(torch.tensor(-time_step/tau_mem)))\n",
        "\n",
        "# Number of time steps to simulate\n",
        "num_steps = 200\n",
        "\n",
        "print(f\"Membrane decay rate ('beta'): {beta}\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qi5EnND98cz3"
      },
      "source": [
        "Instantiating Lapicque's neuron only requires the following line of code:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wI2Nahsg8d-t"
      },
      "source": [
        "# leaky integrate and fire neuron\n",
        "lif1 = snn.Lapicque(beta=beta)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AUv8_QPCIsDT"
      },
      "source": [
        "The same thing can also be accomplished by specifying the RC values:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Q0tk_CBoIv6d"
      },
      "source": [
        "R = 5\n",
        "C = 1e-3\n",
        "\n",
        "lif1 = snn.Lapicque(R=R, C=C, time_step=time_step)\n",
        "\n",
        "print(f\"Membrane decay rate ('beta'): {lif1.beta[0]}\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kuytxXet8lc9"
      },
      "source": [
        "To use this neuron: \n",
        "\n",
        "**Inputs**\n",
        "* `spk_in`: each element of $I_{\\rm in}$, which are all `0` for now, is sequentially passed as an input\n",
        "* `mem`: the membrane potential at the present time $t$ is also passed as input. Initialize it arbitrarily as $U_0 = 0.9~V$.\n",
        "\n",
        "**Outputs**\n",
        "* `spk_out`: output spike $S_{\\rm out}[t+1]$ at the next time step ('1' if there is a spike; '0' if there is no spike)\n",
        "* `mem`: membrane potential $U_{\\rm mem}[t+1]$ at the next time step\n",
        "\n",
        "These all need to be of type `torch.Tensor`.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "M_0O1Q3Y9KgM"
      },
      "source": [
        "# Initialize membrane, input, and output\n",
        "mem = torch.ones(1) * 0.9  # membrane potential of 0.9 at t=0\n",
        "cur_in = torch.zeros(num_steps)  # input is 0 for all t \n",
        "spk_out = torch.zeros(1)  # neuron needs somewhere to sequentially dump its output spikes"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_CEuNHsN8-67"
      },
      "source": [
        "These values are only for the initial time step $t=0$. We'd like to watch the evolution of `mem` over time. The list `mem_rec` is initialized to record these values at every time step."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CUCxQBTQ9d_P"
      },
      "source": [
        "# Initialize somewhere to store recordings of membrane potential\n",
        "mem_rec = [mem]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i4TxVFaV9uf6"
      },
      "source": [
        "Now it's time to run a simulation! 200 time steps will be simulated, updating `mem` at each step and recording its value in `mem_rec`:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "j6GeMMiU8kc4"
      },
      "source": [
        "# pass updated value of mem and cur_in[step]=0 at every time step\n",
        "for step in range(num_steps):\n",
        "  spk_out, mem = lif1(cur_in[step], mem)\n",
        "\n",
        "  # Store recordings of membrane potential\n",
        "  mem_rec.append(mem)\n",
        "\n",
        "# crunch the list of tensors into one tensor\n",
        "mem_rec = torch.stack(mem_rec)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4lngTsykIc7k"
      },
      "source": [
        "Let's take a look at how the membrane potential and synaptic current evolved."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SrfeQWT6I2JC"
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "\n",
        "plt.title(\"Lapicque's Neuron Model Without Stimulus\")\n",
        "plt.plot(mem_rec, label=\"Membrane Potential\")\n",
        "plt.xlabel(\"Time step\")\n",
        "plt.ylabel(\"Membrane Potential\")\n",
        "plt.xlim([0, 50])\n",
        "plt.ylim([0, 1])\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9k0PKng2JYcx"
      },
      "source": [
        "This matches the dynamics that were previously derived. We've shown ourselves that the membrane potential will decay over time in the absence of any input stimuli. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7Z01CioKJpkr"
      },
      "source": [
        "### 2.3.2 Lapicque: Step Input"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B6GZUEFgJ0l0"
      },
      "source": [
        "Now let's apply a step current $I_{\\rm in}(t)$ that switches on at $t=t_0$. Given the linear first-order differential equation:\n",
        "\n",
        "$$ \\tau \\frac{dU_{\\rm mem}}{dt} = -U_{\\rm mem} + RI_{\\rm in}(t),$$\n",
        "\n",
        "the general solution will be:\n",
        "\n",
        "$$U_{\\rm mem}=I_{\\rm in}(t)R + [U_0 - I_{\\rm in}(t)R]e^{-\\frac{t}{\\tau}}$$\n",
        "\n",
        "If the membrane potential is initialized to $U_{\\rm mem}(t=0) = 0 V$, then:\n",
        "\n",
        "$$U_{\\rm mem}(t)=I_{\\rm in}(t)R [1 - e^{-\\frac{t}{\\tau}}]$$\n",
        "\n",
        "Let's visualize what this looks like by triggering a current pulse of $I_{in}=100mA$ at $t_0 = 10ms$."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kD7cLKSuLu9n"
      },
      "source": [
        "# Initialize input current pulse\n",
        "cur_in = torch.cat((torch.zeros(10), torch.ones(190)*0.1), 0)  # input current turns on at t=10\n",
        "\n",
        "# Initialize membrane and output\n",
        "mem = torch.zeros(1)  # membrane potential of 0 at t=0\n",
        "spk_out = torch.zeros(1)  # neuron needs somewhere to sequentially dump its output spikes\n",
        "\n",
        "# Initialize somewhere to store recordings of membrane potential\n",
        "mem_rec = [mem]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PW8ZDbKCNW8E"
      },
      "source": [
        "As before, 200 time steps will be simulated. But this time, the new values of `cur_in` will be passed:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "w9J2AsrZNXr8"
      },
      "source": [
        "# pass updated value of mem and cur_in[step] at every time step\n",
        "for step in range(num_steps):\n",
        "  spk_out, mem = lif1(cur_in[step], mem)\n",
        "\n",
        "  # Store recordings of membrane potential\n",
        "  mem_rec.append(mem)\n",
        "\n",
        "# crunch -list- of tensors into one tensor\n",
        "mem_rec = torch.stack(mem_rec)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_cC3zoy7OW_Y"
      },
      "source": [
        "fig, ax = plt.subplots(2, figsize=(8,6),sharex=True)\n",
        "\n",
        "# Plot input current\n",
        "ax[0].plot(cur_in, c=\"tab:orange\")\n",
        "ax[0].set_ylim([0, 0.2])\n",
        "ax[0].set_ylabel(\"Input Current ($I_{in}$)\")\n",
        "ax[0].set_title(\"Lapicque's Neuron Model With Step Input\")\n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(mem_rec)\n",
        "ax[1].set_ylim([0, 0.6])\n",
        "ax[1].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "\n",
        "ax[1].axvline(x=10, ymin=0, ymax=2.2, alpha = 0.25, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uZi44ZEvNhkR"
      },
      "source": [
        "The membrane potential exponentially rises and then stabilizes at $U_{\\rm mem}(t\\rightarrow \\infty) =I_{\\rm in}R$:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1l40Efnbk92o"
      },
      "source": [
        "print(f\"The calculated value of input pulse [A] x resistance [Ω] is: {cur_in[11]*lif1.R} V\")\n",
        "print(f\"The simulated value of steady-state membrane potential is: {mem_rec[200][0]} V\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "H0ql9cAzpN5D"
      },
      "source": [
        "Close enough!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZqEtGdKAaIbK"
      },
      "source": [
        "### 2.3.3 Lapicque: Pulse Input"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HOgCeLPYaLKZ"
      },
      "source": [
        "Now what if the step input was clipped at $t=30ms$?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hjtupQNxaQ5D"
      },
      "source": [
        "# Initialize input current pulse.\n",
        "cur_in1 = torch.cat((torch.zeros(10), torch.ones(20)*(0.1), torch.zeros(170)), 0)  # input turns on at t=10, off at t=30\n",
        "\n",
        "# Initialize membrane and output\n",
        "mem = torch.zeros(1)  # membrane potential of 0 at t=0\n",
        "spk_out = torch.zeros(1)  # neuron needs somewhere to sequentially dump its output spikes\n",
        "\n",
        "# Initialize somewhere to store recordings of membrane potential\n",
        "mem_rec1 = [mem]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "O8JFL9CGa0NW"
      },
      "source": [
        "# pass updated value of mem and cur_in[step] at every time step\n",
        "for step in range(num_steps):\n",
        "  spk_out, mem = lif1(cur_in1[step], mem)\n",
        "\n",
        "  # Store recordings of membrane potential\n",
        "  mem_rec1.append(mem)\n",
        "\n",
        "# crunch -list- of tensors into one tensor\n",
        "mem_rec1 = torch.stack(mem_rec1)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UM8GGzWOa3GE"
      },
      "source": [
        "fig, ax = plt.subplots(2, figsize=(8,6),sharex=True)\n",
        "\n",
        "# Plot input current\n",
        "ax[0].plot(cur_in1, c=\"tab:orange\")\n",
        "ax[0].set_ylim([0, 0.2])\n",
        "ax[0].set_ylabel(\"Input Current ($I_{in}$)\")\n",
        "ax[0].set_title(\"Lapicque's Neuron Model With Input Pulse\")\n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(mem_rec1)\n",
        "ax[1].set_ylim([0, 1])\n",
        "ax[1].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "\n",
        "ax[1].axvline(x=10, ymin=0, ymax=2.2, alpha = 0.25, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "ax[1].axvline(x=30, ymin=0, ymax=2.2, alpha = 0.25, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "orJvjev6a9tR"
      },
      "source": [
        "It appears to rise just as it did for the step input, but now it decays with a time constant of $\\tau$ as in our first simulation. \n",
        "\n",
        "Let's deliver approximately the same amount of charge $Q = I \\times t$ to the circuit in half the time. This means our input current amplitude will need to be increased by a little, and the time window will be decreased."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "j_G_gWTzbo4J"
      },
      "source": [
        "# Increase amplitude of current pulse; half the time.\n",
        "cur_in2 = torch.cat((torch.zeros(10), torch.ones(10)*0.111, torch.zeros(180)), 0)  # input turns on at t=10, off at t=20\n",
        "\n",
        "# Initialize membrane and output\n",
        "mem = torch.zeros(1)  # membrane potential of 0 at t=0\n",
        "spk_out = torch.zeros(1)  # neuron needs somewhere to sequentially dump its output spikes\n",
        "\n",
        "# Initialize somewhere to store recordings of membrane potential\n",
        "mem_rec2 = [mem]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4YYr91yzbq-G"
      },
      "source": [
        "# pass updated value of mem and cur_in[step] at every time step\n",
        "for step in range(num_steps):\n",
        "  spk_out, mem = lif1(cur_in2[step], mem)\n",
        "\n",
        "  # Store recordings of membrane potential\n",
        "  mem_rec2.append(mem)\n",
        "\n",
        "# crunch -list- of tensors into one tensor\n",
        "mem_rec2 = torch.stack(mem_rec2)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "krBTZm92btXH"
      },
      "source": [
        "fig, ax = plt.subplots(2, figsize=(8,6),sharex=True)\n",
        "\n",
        "# Plot input current\n",
        "ax[0].plot(cur_in2, c=\"tab:orange\")\n",
        "ax[0].set_ylim([0, 0.2])\n",
        "ax[0].set_ylabel(\"Input Current ($I_{in}$)\")\n",
        "ax[0].set_title(\"Lapicque's Neuron Model With Input Pulse: x1/2 pulse width\")\n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(mem_rec2)\n",
        "ax[1].set_ylim([0, 1])\n",
        "ax[1].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "\n",
        "ax[1].axvline(x=10, ymin=0, ymax=2.2, alpha = 0.25, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "ax[1].axvline(x=20, ymin=0, ymax=2.2, alpha = 0.25, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "34azxvjWbqge"
      },
      "source": [
        "Let's do that again, but with an even faster input pulse and higher amplitude:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "650G94v9dgr3"
      },
      "source": [
        "# Increase amplitude of current pulse; quarter the time.\n",
        "cur_in3 = torch.cat((torch.zeros(10), torch.ones(5)*0.147, torch.zeros(185)), 0)  # input turns on at t=10, off at t=15\n",
        "\n",
        "# Initialize membrane and output\n",
        "mem = torch.zeros(1)  # membrane potential of 0 at t=0\n",
        "spk_out = torch.zeros(1)  # neuron needs somewhere to sequentially dump its output spikes\n",
        "\n",
        "# Initialize somewhere to store recordings of membrane potential\n",
        "mem_rec3 = [mem]\n",
        "\n",
        "# pass updated value of mem and cur_in[step] at every time step\n",
        "for step in range(num_steps):\n",
        "  spk_out, mem = lif1(cur_in3[step], mem)\n",
        "\n",
        "  # Store recordings of membrane potential\n",
        "  mem_rec3.append(mem)\n",
        "\n",
        "# crunch -list- of tensors into one tensor\n",
        "mem_rec3 = torch.stack(mem_rec3)\n",
        "\n",
        "# Generate Plots\n",
        "fig, ax = plt.subplots(2, figsize=(8,6),sharex=True)\n",
        "\n",
        "# Plot input current\n",
        "ax[0].plot(cur_in3, c=\"tab:orange\")\n",
        "ax[0].set_ylim([0, 0.2])\n",
        "ax[0].set_ylabel(\"Input Current ($I_{in}$)\")\n",
        "ax[0].set_title(\"Lapicque's Neuron Model With Input Pulse: x1/4 pulse width\")\n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(mem_rec3)\n",
        "ax[1].set_ylim([0, 1])\n",
        "ax[1].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "\n",
        "ax[1].axvline(x=10, ymin=0, ymax=2.2, alpha = 0.25, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "ax[1].axvline(x=15, ymin=0, ymax=2.2, alpha = 0.25, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fuQtNgw4sYhk"
      },
      "source": [
        "Let's compare all three experiments on the same plot:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZL7P3mEfsdAR"
      },
      "source": [
        "# Generate Plots\n",
        "fig, ax = plt.subplots(2, figsize=(8,6),sharex=True)\n",
        "\n",
        "# Plot input current\n",
        "ax[0].plot(cur_in1)\n",
        "ax[0].plot(cur_in2)\n",
        "ax[0].plot(cur_in3)\n",
        "ax[0].set_ylim([0, 0.2])\n",
        "ax[0].set_ylabel(\"Input Current ($I_{in}$)\")\n",
        "ax[0].set_title(\"Lapicque's Neuron Model With Input Pulse: Varying inputs\")\n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(mem_rec1)\n",
        "ax[1].plot(mem_rec2)\n",
        "ax[1].plot(mem_rec3)\n",
        "ax[1].set_ylim([0, 1])\n",
        "ax[1].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "\n",
        "ax[1].axvline(x=10, ymin=0, ymax=2.2, alpha = 0.25, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "ax[1].axvline(x=15, ymin=0, ymax=2.2, alpha = 0.25, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "ax[1].axvline(x=20, ymin=0, ymax=2.2, alpha = 0.25, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "ax[1].axvline(x=30, ymin=0, ymax=2.2, alpha = 0.25, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6WpZx0NStNdd"
      },
      "source": [
        "As the input current pulse amplitude increases, the rise time of the membrane potential speeds up. In the limit of the input current pulse width becoming infinitesimally small, $T_W \\rightarrow 0s$, the membrane potential will jump straight up in virtually zero rise time:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "h0EG17NpuzuN"
      },
      "source": [
        "# Current spike input\n",
        "cur_in4 = torch.cat((torch.zeros(10), torch.ones(1)*0.5, torch.zeros(189)), 0)  # input only on for 1 time step\n",
        "\n",
        "# Initialize membrane and output\n",
        "mem = torch.zeros(1)  # membrane potential of 0 at t=0\n",
        "spk_out = torch.zeros(1)  # neuron needs somewhere to sequentially dump its output spikes\n",
        "\n",
        "# Initialize somewhere to store recordings of membrane potential\n",
        "mem_rec4 = [mem]\n",
        "\n",
        "# pass updated value of mem and cur_in[step] at every time step\n",
        "for step in range(num_steps):\n",
        "  spk_out, mem = lif1(cur_in4[step], mem)\n",
        "\n",
        "  # Store recordings of membrane potential\n",
        "  mem_rec4.append(mem)\n",
        "\n",
        "# crunch -list- of tensors into one tensor\n",
        "mem_rec4 = torch.stack(mem_rec4)\n",
        "\n",
        "# Generate Plots\n",
        "fig, ax = plt.subplots(2, figsize=(8,6),sharex=True)\n",
        "\n",
        "# Plot input current\n",
        "ax[0].plot(cur_in4, c=\"tab:orange\")\n",
        "ax[0].set_ylim([0, 0.6])\n",
        "ax[0].set_ylabel(\"Input Current ($I_{in}$)\")\n",
        "ax[0].set_title(\"Lapicque's Neuron Model With Input Spike\")\n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(mem_rec4)\n",
        "ax[1].set_ylim([0, 1])\n",
        "ax[1].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "ax[1].axvline(x=10, ymin=0, ymax=2.2, alpha = 0.25, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DNY_EftJvREj"
      },
      "source": [
        "The current pulse width is now so short, it effectively looks like a spike. That is to say, charge is delivered in an infinitely short period of time, $I_{\\rm in}(t) = Q/t_0$ where $t_0 \\rightarrow 0$. More formally:\n",
        "\n",
        "$$I_{\\rm in}(t) = Q \\delta (t-t_0),$$\n",
        "\n",
        "where $\\delta (t-t_0)$ is the Dirac-Delta function. Physically, it is impossible to 'instantaneously' deposit charge. But integrating $I_{\\rm in}$ gives a result that makes physical sense, as we can obtain the charge delivered:\n",
        "\n",
        "$$1 = \\int^{t_0 + a}_{t_0 - a}\\delta(t-t_0)dt$$\n",
        "\n",
        "$$f(t_0) = \\int^{t_0 + a}_{t_0 - a}f(t)\\delta(t-t_0)dt$$\n",
        "\n",
        "Here, $f(t_0) = I_{\\rm in}(t_0=10) = 0.5A \\implies f(t) = Q = 0.5C$.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Danwy1q9xfvr"
      },
      "source": [
        "### 2.3.4 Lapicque: Firing"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CUH8IYh9xlV-"
      },
      "source": [
        "So far, we have only seen how a neuron will react to spikes at the input. For a neuron to generate and emit its own spikes at the output, we need to combine the passive membrane model with a threshold.\n",
        "\n",
        "If the membrane potential exceeds this threshold, then a voltage spike will be generated, external to the passive membrane model. \n",
        "\n",
        "<center>\n",
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/examples/tutorial2/2_4_spiking.png?raw=true' width=\"450\">\n",
        "</center>\n",
        "\n",
        "By default, `threshold=1` for all neuron models in snnTorch. So applying a step current input that is insufficient will result in the neuron to function only in the subthreshold regime. This time, we will create a list called `spk_rec` to record any output spikes if they occur. The current step will be set to $I_{\\rm in} = 0.15 A$. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Fe61PZ33w49_"
      },
      "source": [
        "# Small step current input\n",
        "cur_in = torch.cat((torch.zeros(10), torch.ones(190)*0.15), 0)\n",
        "\n",
        "# Initialize membrane and output\n",
        "mem = torch.zeros(1)\n",
        "spk_out = torch.zeros(1) \n",
        "mem_rec = [mem]\n",
        "spk_rec = [spk_out]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HPx-uuH90JTo"
      },
      "source": [
        "# Create a new neuron with a slow time constant\n",
        "lif2 = snn.Lapicque(R=5, C=10)\n",
        "\n",
        "print(f\"Membrane decay rate ('beta'): {lif2.beta[0]}\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "se0fnQ5EDHr9"
      },
      "source": [
        "Note how this new value of $\\beta$ is much larger than `lif1.beta`$\\approx 0.82$. \n",
        "\n",
        "For `lif2.beta`$=0.98$, the membrane potential is 98% of the value of that of the previous time step, and experiences a much slower decay rate."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0yqODQazwg4s"
      },
      "source": [
        "# Simulation run across 200 time steps. \n",
        "for step in range(num_steps):\n",
        "  spk_out, mem = lif2(cur_in[step], mem)\n",
        "\n",
        "  # record outputs over time\n",
        "  mem_rec.append(mem)\n",
        "  spk_rec.append(spk_out)\n",
        "\n",
        "# crunch -list- of tensors into one tensor\n",
        "mem_rec = torch.stack(mem_rec)\n",
        "spk_rec = torch.stack(spk_rec)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Z3Laa83DxPbJ"
      },
      "source": [
        "# Generate Plots\n",
        "fig, ax = plt.subplots(2, figsize=(8,6),sharex=True)\n",
        "\n",
        "# Plot input current\n",
        "ax[0].plot(cur_in, c=\"tab:orange\")\n",
        "ax[0].set_ylim([0, 0.4])\n",
        "ax[0].set_ylabel(\"Input Current ($I_{in}$)\")\n",
        "ax[0].set_title(\"Lapicque's Neuron Model With Step Input\")\n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(mem_rec)\n",
        "ax[1].set_ylim([0, 1.25])\n",
        "ax[1].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "ax[1].axhline(y=1.0, alpha=0.25, linestyle=\"dashed\", c=\"black\", linewidth=2)\n",
        "\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vrAD1fBextPM"
      },
      "source": [
        "The membrane potential fails to reach the threshold of 1.0. Instead, it reaches the steady-state value of $I_{\\rm in}R = 0.15A \\times 5Ω = 0.75V$:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XtkMfF5Qxrvm"
      },
      "source": [
        "print(f\"The calculated steady state membrane potential is: {lif1.R*cur_in[199]}\")\n",
        "print(f\"The simulated steady state membrane potential is: {mem_rec[199][0]}\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FN1Gmvuwy1ZA"
      },
      "source": [
        "> Note: these are non-biologically accurate values, and are chosen for simplicity.\n",
        "\n",
        "To reach the threshold, we need to ensure that $I_{\\rm in}R > U_{\\rm thr}$. So set $I_{\\rm in} = 0.21 A$:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LU9dSGi6zEOP"
      },
      "source": [
        "# Larger current step\n",
        "cur_in = torch.cat((torch.zeros(10), torch.ones(190)*0.21), 0)\n",
        "\n",
        "# Initialize membrane and output\n",
        "mem = torch.zeros(1)\n",
        "spk_out = torch.zeros(1) \n",
        "mem_rec = [mem]\n",
        "spk_rec = [spk_out]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "td4oIyzNzLft"
      },
      "source": [
        "# Simulation run across 200 time steps.\n",
        "for step in range(num_steps):\n",
        "  spk_out, mem = lif2(cur_in[step], mem)\n",
        "\n",
        "  # record outputs over time\n",
        "  mem_rec.append(mem)\n",
        "  spk_rec.append(spk_out)\n",
        "\n",
        "# crunch -list- of tensors into one tensor\n",
        "mem_rec = torch.stack(mem_rec)\n",
        "spk_rec = torch.stack(spk_rec)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IeiswwDXGnay"
      },
      "source": [
        "To plot our results, let's import `snntorch.spikeplot`. From Tutorial 1, we learnt how to use it to create raster plots of spike responses."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gxLfSFTIGvaU"
      },
      "source": [
        "from snntorch import spikeplot as splt"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4BYXOPSizWH6"
      },
      "source": [
        "# Generate Plots\n",
        "fig, ax = plt.subplots(3, figsize=(8,6), sharex=True, \n",
        "                       gridspec_kw = {'height_ratios': [1, 1, 0.4]})\n",
        "\n",
        "# Plot input current\n",
        "ax[0].plot(cur_in, c=\"tab:orange\")\n",
        "ax[0].set_ylim([0, 0.4])\n",
        "ax[0].set_xlim([0, 200])\n",
        "ax[0].set_ylabel(\"Input Current ($I_{in}$)\")\n",
        "ax[0].set_title(\"Lapicque's Neuron Model With Step Input\")\n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(mem_rec)\n",
        "ax[1].set_ylim([0, 1.25])\n",
        "ax[1].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "ax[1].axhline(y=1.0, alpha=0.25, linestyle=\"dashed\", c=\"black\", linewidth=2)\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "# Plot output spike using spikeplot\n",
        "splt.raster(spk_rec, ax[2], s=400, c=\"black\", marker=\"|\")\n",
        "ax[2].axvline(x=162, ymin=0, ymax=6.75, alpha = 0.15, linestyle=\"dashed\", c=\"black\", linewidth=2, zorder=0, clip_on=False)\n",
        "plt.ylabel(\"Output spikes\")\n",
        "plt.yticks([]) \n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cm5tEp4YEoQm"
      },
      "source": [
        "The membrane potential exponentially rises and then hits the threshold, at which point it resets. We can roughly see this occurs between $155s < t_{\\rm spk} < 165s$:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tWSGUUsgNr9D"
      },
      "source": [
        "print(spk_rec[155:165])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4xWW-ncoOWrJ"
      },
      "source": [
        "The absence of a spike is represented by $S_{\\rm out}=0$, and the occurrence of a spike is $S_{\\rm out}=1$. Here, the spike occurs at $S_{\\rm out}(t=162)=1$.\n",
        "\n",
        "If you are wondering why each of these entries is stored as a tensor, it is because soon we will simulate large scale neural networks. Each entry will contain the spike responses of many neurons, and tensors can be loaded into GPU memory to speed up the training process.\n",
        "\n",
        "If $I_{\\rm in}$ is increased, then the membrane potential approaches $U_{\\rm thr}$ faster:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HY6G9HZ1POC7"
      },
      "source": [
        "# Even Larger current step\n",
        "cur_in = torch.cat((torch.zeros(10), torch.ones(190)*0.3), 0)\n",
        "\n",
        "# Initialize membrane and output\n",
        "mem = torch.zeros(1)\n",
        "spk_out = torch.zeros(1) \n",
        "mem_rec = [mem]\n",
        "spk_rec = [spk_out]\n",
        "\n",
        "# Simulation run across 200 time steps.\n",
        "for step in range(num_steps):\n",
        "  spk_out, mem = lif2(cur_in[step], mem)\n",
        "\n",
        "  # record outputs over time\n",
        "  mem_rec.append(mem)\n",
        "  spk_rec.append(spk_out)\n",
        "\n",
        "# crunch -list- of tensors into one tensor\n",
        "mem_rec = torch.stack(mem_rec)\n",
        "spk_rec = torch.stack(spk_rec)\n",
        "\n",
        "# Generate Plots\n",
        "fig, ax = plt.subplots(3, figsize=(8,6), sharex=True, \n",
        "                       gridspec_kw = {'height_ratios': [1, 1, 0.4]})\n",
        "\n",
        "# Plot input current\n",
        "ax[0].plot(cur_in, c=\"tab:orange\")\n",
        "ax[0].set_ylim([0, 0.4])\n",
        "ax[0].set_xlim([0, 200])\n",
        "ax[0].set_ylabel(\"Input Current ($I_{in}$)\")\n",
        "ax[0].set_title(\"Lapicque's Neuron Model With Periodic Firing\")\n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(mem_rec)\n",
        "ax[1].set_ylim([0, 1.25])\n",
        "ax[1].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "ax[1].axhline(y=1.0, alpha=0.25, linestyle=\"dashed\", c=\"black\", linewidth=2)\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "# Plot output spike using spikeplot\n",
        "splt.raster(spk_rec, ax[2], s=400, c=\"black\", marker=\"|\")\n",
        "plt.ylabel(\"Output spikes\")\n",
        "plt.yticks([]) \n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GOjWXepnPleo"
      },
      "source": [
        "A similar increase in firing frequency can also be induced by decreasing the threshold. This requires initializing a new neuron model, but the rest of the code block is the exact same as above:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HQ0VIBlJPrvr"
      },
      "source": [
        "# Half the threshold\n",
        "lif3 = snn.Lapicque(R=5, C=10, threshold=0.5)\n",
        "\n",
        "# Initialize membrane and output\n",
        "mem = torch.zeros(1)\n",
        "spk_out = torch.zeros(1) \n",
        "mem_rec = [mem]\n",
        "spk_rec = [spk_out]\n",
        "\n",
        "# Simulation run across 200 time steps.\n",
        "for step in range(num_steps):\n",
        "  spk_out, mem = lif3(cur_in[step], mem)\n",
        "\n",
        "  # record outputs over time\n",
        "  mem_rec.append(mem)\n",
        "  spk_rec.append(spk_out)\n",
        "\n",
        "# crunch -list- of tensors into one tensor\n",
        "mem_rec = torch.stack(mem_rec)\n",
        "spk_rec = torch.stack(spk_rec)\n",
        "\n",
        "# Generate Plots\n",
        "fig, ax = plt.subplots(3, figsize=(8,6), sharex=True, \n",
        "                       gridspec_kw = {'height_ratios': [1, 1, 0.4]})\n",
        "\n",
        "# Plot input current\n",
        "ax[0].plot(cur_in, c=\"tab:orange\")\n",
        "ax[0].set_ylim([0, 0.4])\n",
        "ax[0].set_xlim([0, 200])\n",
        "ax[0].set_ylabel(\"Input Current ($I_{in}$)\")\n",
        "ax[0].set_title(\"Lapicque's Neuron Model With Low Threshold\")\n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(mem_rec)\n",
        "ax[1].set_ylim([0, 1.25])\n",
        "ax[1].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "ax[1].axhline(y=0.5, alpha=0.25, linestyle=\"dashed\", c=\"black\", linewidth=2)\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "# Plot output spike using spikeplot\n",
        "splt.raster(spk_rec, ax[2], s=400, c=\"black\", marker=\"|\")\n",
        "plt.ylabel(\"Output spikes\")\n",
        "plt.yticks([]) \n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rJwtafQ_Siav"
      },
      "source": [
        "That's what happens for a constant current injection. But in both deep neural networks and in the biological brain, most neurons will be connected to other neurons. They are more likely to receive spikes, rather than injections of constant current. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sOr9akj9S2Ei"
      },
      "source": [
        "### 2.3.5 Lapicque: Spike Inputs"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rw_FZDnMS-pc"
      },
      "source": [
        "Let's harness some of the skills we learnt in [Tutorial 1](https://colab.research.google.com/github/jeshraghian/snntorch/blob/tutorials/examples/tutorial_1_spikegen.ipynb), and use the `snntorch.spikegen` module to create some randomly generated input spikes."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "V7AXG4QqTCoZ"
      },
      "source": [
        "from snntorch import spikegen \n",
        "\n",
        "# Create a 1-D random spike train. Each element has a probability of 40% of firing.\n",
        "spk_in = spikegen.rate_conv(torch.ones((num_steps)) * 0.40)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5CJxZg1TTNCJ"
      },
      "source": [
        "Run the following code block to see how many spikes have been generated."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rbhLmpNSTPgK"
      },
      "source": [
        "# Tell me the number of spikes\n",
        "print(f\"There are {int(sum(spk_in))} total spikes out of {len(spk_in)} time steps.\")\n",
        "\n",
        "# Now show me the spikes\n",
        "from snntorch import spikeplot as splt\n",
        "\n",
        "fig = plt.figure(facecolor=\"w\", figsize=(8, 1))\n",
        "ax = fig.add_subplot(111)\n",
        "\n",
        "splt.raster(spk_in.reshape(num_steps, -1), ax, s=100, c=\"black\", marker=\"|\")\n",
        "\n",
        "plt.title(\"Input Spikes\")\n",
        "plt.xlabel(\"Time step\")\n",
        "plt.yticks([])\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wo2_Voc3TvSc"
      },
      "source": [
        "# Refresh all our hidden and output variables\n",
        "mem = torch.ones(1)*0.5  # membrane potential of 0.5 at t=0\n",
        "spk_out = torch.zeros(1)  # neuron needs somewhere to dump its output spikes\n",
        "\n",
        "# Create a trace of the variables of interest\n",
        "mem_rec = [mem]\n",
        "spk_rec = [spk_out]\n",
        "\n",
        "# Run the simulation\n",
        "for step in range(num_steps):\n",
        "  spk_out, mem = lif3(spk_in[step], mem)\n",
        "\n",
        "  # Store recordings of output and hidden states\n",
        "  spk_rec.append(spk_out)\n",
        "  mem_rec.append(mem)\n",
        "\n",
        "# crunch -list- of tensors into one tensor\n",
        "mem_rec = torch.stack(mem_rec)\n",
        "spk_rec = torch.stack(spk_rec)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uGYz-PJAT_1q"
      },
      "source": [
        "# Generate Plots\n",
        "fig, ax = plt.subplots(3, figsize=(8,6), sharex=True, \n",
        "                       gridspec_kw = {'height_ratios': [0.4, 1, 0.4]})\n",
        "\n",
        "# Plot input current\n",
        "splt.raster(spk_in, ax[0], s=400, c=\"black\", marker=\"|\")\n",
        "ax[0].set_ylabel(\"Input Spikes\")\n",
        "ax[0].set_title(\"Lapicque's Neuron Model With Input Spikes\")\n",
        "plt.yticks([]) \n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(mem_rec)\n",
        "ax[1].set_ylim([0, 1])\n",
        "ax[1].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "ax[1].axhline(y=0.5, alpha=0.25, linestyle=\"dashed\", c=\"black\", linewidth=2)\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "# Plot output spike using spikeplot\n",
        "splt.raster(spk_rec, ax[2], s=400, c=\"black\", marker=\"|\")\n",
        "plt.ylabel(\"Output spikes\")\n",
        "plt.yticks([]) \n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fGYnDIoTXBCo"
      },
      "source": [
        "### 2.3.6 Lapicque: Reset Mechanisms"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lQsH8VOGXF2Q"
      },
      "source": [
        "The final detail of the Lapicque neuron we want to explore is the sharp drop of membrane potential every time the neuron emits an output spike. This sharp drops promotes a reduction of spike generation, which supplements part of the theory on how brains are so power efficient. Biologically, this is known as the 'refractory period' where the the neuron's firing ability is momentarily suppressed. Here, we use a reset mechanism to model the refractory period.\n",
        "\n",
        "There are two ways to implement the reset mechanism:\n",
        "\n",
        "1.  *reset by subtraction* (default) $-$ subtract the threshold from the membrane potential each time a spike is generated;\n",
        "2.   *reset to zero* $-$ force the membrane potential to zero each time a spike is generated.\n",
        "\n",
        "<center>\n",
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/examples/tutorial2/2_5_reset.png?raw=true' width=\"450\">\n",
        "</center>\n",
        "\n",
        "Let's instantiate another neuron model to demonstrate how to alternate between reset mechanisms. \n",
        "\n",
        "By default, snnTorch neuron models use `reset_mechanism = \"subtract\"`. This can be explicitly overridden by passing the argument `reset_mechanism =  \"zero\"`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7bfidqhKXdri"
      },
      "source": [
        "# Reset mechanism: zero\n",
        "lif4 = snn.Lapicque(R=5, C=10, threshold=0.5, reset_mechanism=\"zero\")\n",
        "\n",
        "# Refresh all our hidden and output variables\n",
        "mem = torch.ones(1)*0.5  # membrane potential of 0.5 at t=0\n",
        "spk_out = torch.zeros(1)  # neuron needs somewhere to dump its output spikes\n",
        "\n",
        "# Create a trace of the variables of interest\n",
        "mem_rec0 = [mem]\n",
        "spk_rec0 = [spk_out]\n",
        "\n",
        "# Run the simulation\n",
        "for step in range(num_steps):\n",
        "  spk_out, mem = lif4(spk_in[step], mem)\n",
        "\n",
        "  # Store recordings of output and hidden states\n",
        "  spk_rec0.append(spk_out)\n",
        "  mem_rec0.append(mem)\n",
        "\n",
        "# crunch -list- of tensors into one tensor\n",
        "mem_rec0 = torch.stack(mem_rec0)\n",
        "spk_rec0 = torch.stack(spk_rec0)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ifU3betRYETs"
      },
      "source": [
        "# Generate Plots to Compare Reset Mechanisms\n",
        "fig, ax = plt.subplots(nrows=3, ncols=2, figsize=(10,6), sharex=True, \n",
        "                       gridspec_kw = {'height_ratios': [0.4, 1, 0.4], 'wspace':0.05})\n",
        "\n",
        "# Reset by Subtraction: input spikes\n",
        "splt.raster(spk_in, ax[0][0], s=400, c=\"black\", marker=\"|\")\n",
        "ax[0][0].set_ylabel(\"Input Spikes\")\n",
        "ax[0][0].set_title(\"Reset by Subtraction\")\n",
        "ax[0][0].set_yticks([])\n",
        "\n",
        "# Reset by Subtraction: membrane potential \n",
        "ax[1][0].plot(mem_rec)\n",
        "ax[1][0].set_ylim([0, 0.7])\n",
        "ax[1][0].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "ax[1][0].axhline(y=0.5, alpha=0.25, linestyle=\"dashed\", c=\"black\", linewidth=2)\n",
        "\n",
        "# Reset by Subtraction: output spikes\n",
        "splt.raster(spk_rec, ax[2][0], s=400, c=\"black\", marker=\"|\")\n",
        "ax[2][0].set_yticks([])\n",
        "ax[2][0].set_xlabel(\"Time step\")\n",
        "ax[2][0].set_ylabel(\"Output Spikes\")\n",
        "\n",
        "# Reset to Zero: input spikes\n",
        "splt.raster(spk_in, ax[0][1], s=400, c=\"black\", marker=\"|\")\n",
        "ax[0][1].set_title(\"Reset to Zero\")\n",
        "ax[0][1].set_yticks([])\n",
        "\n",
        "# Reset to Zero: membrane potential\n",
        "ax[1][1].plot(mem_rec0)\n",
        "ax[1][1].set_ylim([0, 0.7])\n",
        "ax[1][1].axhline(y=0.5, alpha=0.25, linestyle=\"dashed\", c=\"black\", linewidth=2)\n",
        "ax[1][1].set_yticks([])\n",
        "ax[2][1].set_xlabel(\"Time step\")\n",
        "\n",
        "# Reset to Zero: output spikes\n",
        "splt.raster(spk_rec0, ax[2][1], s=400, c=\"black\", marker=\"|\")\n",
        "ax[2][1].set_yticks([])\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Dzs7qX2Kb8j_"
      },
      "source": [
        "Pay close attention to the evolution of the membrane potential, especially in the moments after it reaches the threshold. You may notice that for \"Reset to Zero\", the membrane potential is forced back to zero after each spike.\n",
        "\n",
        "So which one is better? Applying `\"subtract\"` (the default value in `reset_mechanism`) is less lossy, because it does not ignore how much the membrane exceeds the threshold by.\n",
        "\n",
        "On the other hand, applying a hard reset with `\"zero\"` promotes sparsity and potentially less power consumption when running on dedicated neuromorphic hardware. Both options are available for you to experiment with. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3ukHPMnOieVV"
      },
      "source": [
        "## 2.4 Synaptic Conductance-based LIF Neuron Model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9b9GEzOfm_Iz"
      },
      "source": [
        "The passive membrane model allows discrete current spikes to be passed directly into the neuron. In reality, a spike will result in the gradual release of neurotransmitters from the pre-synaptic neuron to the post-synaptic neuron. This model accounts for the gradual temporal dynamics of input current, and is no longer strictly modelling a LIF neuron alone."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YO6I0Wpmo3fx"
      },
      "source": [
        "### 2.4.1 Synaptic Current"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KBuE963K_JmQ"
      },
      "source": [
        "If a pre-synaptic neuron fires, the voltage spike is transmitted down the axon of the neuron. It triggers the vesicles to release neurotransmitters into the synaptic cleft. These activate the post-synaptic receptors, which directly influence the effective current that flows into the post-synaptic neuron. \n",
        "\n",
        "Shown below are two types of excitatory receptors.\n",
        "\n",
        "<center>\n",
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/examples/tutorial2/2_6_synaptic.png?raw=true' width=\"600\">\n",
        "</center>\n",
        "\n",
        "The simplest model of synaptic current assumes an increasing current on a very fast time-scale (or instantaneous), followed by a relatively slow exponential decay. This is very similar to the membrane potential dynamics of Lapicque's model.\n",
        "\n",
        "The synaptic condutance-based neuron model combines the synaptic current dynamics with the passive membrane. It must be instantiated with two input arguments:\n",
        "* $\\alpha$: the decay rate of the synaptic current\n",
        "* $\\beta$: the decay rate of the membrane potential (as with Lapicque)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FjJQIlY9DGtR"
      },
      "source": [
        "# Decay rate of LIF states\n",
        "alpha = 0.9\n",
        "beta = 0.8\n",
        "\n",
        "# Initialize 2nd-order LIF neuron\n",
        "lif5 = snn.Synaptic(alpha=alpha, beta=beta)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ssWDmsKlDc34"
      },
      "source": [
        "Using this neuron is the exact same as Lapcique's neuron, but now with the addition of synaptic current `syn` as an input and output:\n",
        "\n",
        "**Inputs**\n",
        "* `spk_in`: each input voltage spike $S_{\\rm in}[t]$ is sequentially passed in\n",
        "* `syn`: synaptic current $I_{\\rm syn}[t]$ at the present time $t$\n",
        "* `mem`: membrane potential $U_{\\rm mem}[t]$ at the present time $t$\n",
        "\n",
        "**Outputs**\n",
        "* `spk_out`: output spike $S_{\\rm out}[t+1]$ at the next time step ('1' if there is a spike; '0' if there is no spike)\n",
        "* `syn`: synaptic current $I_{\\rm syn}[t+1]$ at the next time step\n",
        "* `mem`: membrane potential $U_{\\rm mem}[t+1]$ at the next time step\n",
        "\n",
        "These all need to be of type `torch.Tensor`."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OCzmKeDRFeBb"
      },
      "source": [
        "Apply a periodic spiking input to see how current and membrane evolve with time:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zWSjrxQGENPE"
      },
      "source": [
        "# Initialize hidden states and output\n",
        "syn = torch.zeros(1)  # synaptic current of 0 at t=0\n",
        "mem = torch.zeros(1)  # membrane potential of 0 at t=0\n",
        "spk_out = torch.zeros(1)  # neuron needs somewhere to dump its output spikes\n",
        "\n",
        "# Periodic spiking input, spk_in = 0.2 V\n",
        "spk_period = torch.cat((torch.ones(1)*0.2, torch.zeros(9)), 0)\n",
        "spk_in = spk_period.repeat(20)\n",
        "\n",
        "syn_rec = [syn]\n",
        "mem_rec = [mem]\n",
        "spk_rec = [spk_out]\n",
        "\n",
        "for step in range(num_steps):\n",
        "  spk_out, syn, mem = lif5(spk_in[step], syn, mem)\n",
        "\n",
        "  # Store recordings of output and hidden states\n",
        "  spk_rec.append(spk_out)\n",
        "  syn_rec.append(syn)\n",
        "  mem_rec.append(mem)\n",
        "\n",
        "# crunch -list- of tensors into one tensor\n",
        "spk_rec = torch.stack(spk_rec)\n",
        "syn_rec = torch.stack(syn_rec)\n",
        "mem_rec = torch.stack(mem_rec)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "b1rVl69UG8xr"
      },
      "source": [
        "# Generate Plots\n",
        "fig, ax = plt.subplots(4, figsize=(8,7), sharex=True, \n",
        "                       gridspec_kw = {'height_ratios': [0.4, 1, 1, 0.4]})\n",
        "\n",
        "# Plot input current\n",
        "splt.raster(spk_in, ax[0], s=400, c=\"black\", marker=\"|\")\n",
        "ax[0].set_ylabel(\"Input Spikes\")\n",
        "ax[0].set_title(\"Synaptic Conductance-based Neuron Model With Input Spikes\")\n",
        "ax[0].set_yticks([]) \n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(syn_rec)\n",
        "ax[1].set_ylim([0, 0.5])\n",
        "ax[1].set_ylabel(\"Synaptic Current ($I_{syn}$)\")\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "# Plot membrane potential\n",
        "ax[2].plot(mem_rec)\n",
        "ax[2].set_ylim([0, 1.5])\n",
        "ax[2].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "ax[2].axhline(y=1, alpha=0.25, linestyle=\"dashed\", c=\"black\", linewidth=2)\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "# Plot output spike using spikeplot\n",
        "splt.raster(spk_rec, ax[3], s=400, c=\"black\", marker=\"|\")\n",
        "plt.ylabel(\"Output spikes\")\n",
        "ax[3].set_yticks([]) \n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_w1uL7kbzfwW"
      },
      "source": [
        "If you're not interested in the mathematical detail, then feel free to skip this brief section. We represent the voltage spikes at the input with:\n",
        "\n",
        "$$S_{\\rm in} = \\sum_k \\delta(t-t_k),$$\n",
        "\n",
        "where each spike triggers a weighted jump in synaptic current at time $t_k$, and is followed by an exponential decay:\n",
        "\n",
        "$$I_{\\rm syn}(t) = \\sum_k W_{i,j} S_{in; i,j}(t) e^{-(t-t_k)/\\tau_{syn}}\\Theta(t-t_k)$$\n",
        "\n",
        "* $W_{i, j}$ is the weight between the the $i^{\\rm th}$ pre-synaptic neuron and the $j^{\\rm th}$ post-synaptic neuron\n",
        "\n",
        "* $t_k$ is the timing of each incident spike\n",
        "\n",
        "* $\\Theta(t)$ is the Heaviside step function, which clips the exponential term such that the contribution from each presynaptic spike commences at $t_k$\n",
        "\n",
        "* $\\tau_{syn}$ is the time constant of the synaptic current, independent of the membrane potential time constant\n",
        "\n",
        "The time constant $\\tau_{syn}$ can be equivalently represented by a decay rate $\\alpha$ that specifies the ratio of synaptic current between subsequent time steps:\n",
        "\n",
        "$$\\alpha = \\frac{e^{-\\frac{1}{\\tau_{syn}}}}{e^{-\\frac{0}{\\tau_{syn}}}} = \\frac{e^{-\\frac{2}{\\tau_{syn}}}}{e^{-\\frac{1}{\\tau_{syn}}}} = \\frac{e^{-\\frac{3}{\\tau_{syn}}}}{e^{-\\frac{2}{\\tau_{syn}}}}=~~...$$\n",
        "$$\\implies \\alpha = e^{-\\frac{1}{\\tau_{syn}}}$$\n",
        "\n",
        "\n",
        "When an input spike arrives at the neuron, the synaptic current will jump up $W_{i,j}S_{\\rm in}(t=t_k)$, where $S_{\\rm in}(t=t_k)=1$. \n",
        "\n",
        "That is to say: $\\Delta I_{\\rm syn}(t=t_k) = W_{i, j}$\n",
        "\n",
        "\n",
        "\n",
        "In summary, each spike contributes a shifted exponential decay to the synaptic current $I_{\\rm syn}$, which are all summed together. This current is then integrated by the passive membrane equation derived in the previous section, thus generating output spikes.\n",
        "\n",
        "If the math doesn't make sense, don't worry about it. A graphical intuition is usually sufficient to understand the essence of the synaptic conductance-based neuron model. \n",
        "\n",
        "<center>\n",
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/examples/tutorial2/2_7_stein.png?raw=true' width=\"600\">\n",
        "</center>\n",
        "\n",
        "This model has the same optional input arguments of `reset_mechanism` and `threshold` as described for Lapicque's neuron model."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "P-m632rfMXIb"
      },
      "source": [
        "# 3. A Feedforward Spiking Neural Network"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iTV4g4JjMgef"
      },
      "source": [
        "So far, we have only considered how one neuron reacts to a single input stimulus. snnTorch makes it extremely straightforward to scale this up to a deep neural network. Here, we will create a 3-layer fully-connected neural network of dimensions 784-1000-10.\n",
        "\n",
        "Compared to our simulations so far, each neuron will now integrate over many more incoming input spikes. \n",
        "\n",
        "\n",
        "<center>\n",
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/examples/tutorial2/2_8_fcn.png?raw=true' width=\"600\">\n",
        "</center>\n",
        "\n",
        "We will use PyTorch to form the connections between neurons, and snnTorch is used to create the neurons."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "whg4LFGFNuH1"
      },
      "source": [
        "import torch \n",
        "import torch.nn as nn"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pypN9pglN2C-"
      },
      "source": [
        "First, initialize all layers."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LDE2ciwCOGa-"
      },
      "source": [
        "num_inputs = 784\n",
        "num_hidden = 1000\n",
        "num_outputs = 10\n",
        "\n",
        "# initialize layers\n",
        "fc1 = nn.Linear(num_inputs, num_hidden)\n",
        "lif1 = snn.Synaptic(alpha=alpha, beta=beta)\n",
        "fc2 = nn.Linear(num_hidden, num_outputs)\n",
        "lif2 = snn.Synaptic(alpha=alpha, beta=beta)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3bh1Op3bOeLk"
      },
      "source": [
        "Next, initialize the hidden variables and outputs of each spiking neuron. \n",
        "As your networks increase in size, this will become a tedious process. So we can call a static method `init_synaptic()` to take care of this. All neurons in snnTorch have their own initialization methods that follow this same syntax, e.g., `init_lapicque()`. The shape of the hidden states are automatically initialized based on the input data dimensions during the first forward pass. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "p3X1_5puO3W7"
      },
      "source": [
        "# Initialize hidden states\n",
        "syn1, mem1 = lif1.init_synaptic()\n",
        "syn2, mem2 = lif2.init_synaptic()\n",
        "\n",
        "# Lists to record output traces\n",
        "mem2_rec = []\n",
        "spk1_rec = []\n",
        "spk2_rec = []"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jotASwAQQmVJ"
      },
      "source": [
        "Create an input spike train to pass into the network. There are 200 time steps to simulate across 784 input neurons. We 'unsqueeze' the input along `dim=1` to denote this to be 'one batch' of data. So the dimensions of this input tensor must be 200 $\\times$ 1 $\\times$ 784:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ON4KKAjTQzf1"
      },
      "source": [
        "spk_in = spikegen.rate_conv(torch.rand((200, 784))*0.1).unsqueeze(1)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CREME2dLPo-G"
      },
      "source": [
        "Now it's finally time to run a full simulation. \n",
        "An intuitive way to think about how PyTorch and snnTorch work together is that PyTorch routes the neurons together, and snnTorch loads the results into spiking neuron models. In terms of coding up a network, these spiking neurons can be treated like time-varying activation functions.\n",
        "\n",
        "Recall that the output of a spiking neuron is $S_{\\rm out}=1$ when a spike is triggered. This spike is then passed to the next layer. It is multiplied by the weight initialized by `nn.Linear` $S_{\\rm out; i}\\times W_{i, j}$, just as the output activation of a standard artificial neuron would be in a non-spiking neural network. The weighted spike is then passed as the input to the next layer of neurons for a given time step. If there is no spike, then nothing is passed to the post-synaptic neuron.\n",
        "\n",
        "The only difference from our simulations thus far is that we sequentially pass the output through additional layers of neurons. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-5Zron7aNkxH"
      },
      "source": [
        "for step in range(num_steps):\n",
        "    cur1 = fc1(spk_in[step])\n",
        "    spk1, syn1, mem1 = lif1(cur1, syn1, mem1)\n",
        "    cur2 = fc2(spk1)\n",
        "    spk2, syn2, mem2 = lif2(cur2, syn2, mem2)\n",
        "\n",
        "    mem2_rec.append(mem2)\n",
        "    spk1_rec.append(spk1)\n",
        "    spk2_rec.append(spk2)\n",
        "\n",
        "# convert output recordings to tensors\n",
        "mem2_rec = torch.stack(mem2_rec)\n",
        "spk1_rec = torch.stack(spk1_rec)\n",
        "spk2_rec = torch.stack(spk2_rec)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xGpWl-S2YAyS"
      },
      "source": [
        "At this stage, the spikes don't have any real meaning. The inputs and weights are all randomly initialized, and no training has taken place. But let's take a look at the raster plots just to check that the spikes are propagating to the output layer."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8FJuQGKfWEmR"
      },
      "source": [
        "# Generate Plots\n",
        "fig, ax = plt.subplots(3, figsize=(8,7), sharex=True, \n",
        "                       gridspec_kw = {'height_ratios': [1, 1, 0.4]})\n",
        "\n",
        "# Plot input spikes\n",
        "splt.raster(spk_in[:,0], ax[0], s=0.05, c=\"black\")\n",
        "ax[0].set_ylabel(\"Input Spikes\")\n",
        "ax[0].set_title(\"Fully Connected Spiking Neural Network\")\n",
        "\n",
        "# Plot hidden layer spikes\n",
        "splt.raster(spk1_rec.reshape(num_steps, -1), ax[1], s = 0.05, c=\"black\")\n",
        "ax[1].set_ylabel(\"Hidden Layer\")\n",
        "\n",
        "# Plot output spikes\n",
        "splt.raster(spk2_rec.reshape(num_steps, -1), ax[2], c=\"black\", marker=\"|\")\n",
        "ax[2].set_ylabel(\"Output Spikes\")\n",
        "ax[2].set_ylim([0, 10])\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jfZwr307P5-u"
      },
      "source": [
        "We can also use `spikeplot.spike_count` to generate a spike counter of the output layer.<br>\n",
        "Note: if you are running the notebook locally on your desktop, please uncomment the line below and modify the path to your ffmpeg.exe"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fAXeoiyGOm1H"
      },
      "source": [
        "from IPython.display import HTML\n",
        "\n",
        "fig, ax = plt.subplots(facecolor='w', figsize=(12, 7))\n",
        "labels=['0', '1', '2', '3', '4', '5', '6', '7', '8','9']\n",
        "spk2_rec = spk2_rec.squeeze(1).detach().cpu()\n",
        "\n",
        "# plt.rcParams['animation.ffmpeg_path'] = 'C:\\\\path\\\\to\\\\your\\\\ffmpeg.exe'\n",
        "\n",
        "#  Plot spike count histogram\n",
        "anim = splt.spike_count(spk2_rec, fig, ax, labels=labels, animate=True)\n",
        "HTML(anim.to_html5_video())\n",
        "# anim.save(\"spike_bar.gif\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xUXpWaPybVkd"
      },
      "source": [
        "We can also visualize the membrane potential traces with `spikeplot.traces`. We'll plot 9 out of 10 output neurons. Compare it to the animation and raster plot above to see if you can match the traces to the neuron. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SpqvImMPxO36"
      },
      "source": [
        "splt.traces(mem2_rec.squeeze(1), spk=spk2_rec.squeeze(1))\n",
        "\n",
        "fig = plt.gcf() \n",
        "fig.set_size_inches(8, 6)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "z98Y9ucNac3e"
      },
      "source": [
        "# 4. Alpha Neuron Model (Hacked Spike Response Model)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "te-QotkqakDL"
      },
      "source": [
        "To finish up this tutorial, a recursive version of the Spike Response Model (SRM), or 'Alpha' neuron, is also available, called using `snntorch.Alpha`. The neuron models thus far have all been based on the passive membrane model, using ordinary differential equations to describe their dynamics.\n",
        "\n",
        "The SRM family of models, on the other hand, is interpreted in terms of a filter. Upon the arrival of an input spike, this spike is convolved with the filter to give the membrane potential response. The form of this filter can be exponential, as is the case with Lapicque's neuron, or they can be more complex such as a sum of exponentials. SRM models are appealing as they can arbitrarily add refractoriness, threshold adaptation, and any number of other features simply by embedding them into the filter. \n",
        "\n",
        "<left>\n",
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/examples/tutorial2/exp.gif?raw=true' width=\"400\">\n",
        "</left> \n",
        "\n",
        "<right>\n",
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/examples/tutorial2/alpha.gif?raw=true' width=\"400\">\n",
        "</right> \n",
        "\n",
        "Formally, this process is represented by:\n",
        "\n",
        "$$U_{\\rm mem}(t) = \\sum_i W_{i, j}(\\epsilon * S_{\\rm in; i,j})(t)$$\n",
        "\n",
        "where the incoming spikes $S_{\\rm in; i,j}$ are convolved with a spike response kernel $\\epsilon( \\cdot )$. The spike response is scaled by a synaptic weight, $W_{i, j}$. In the figures above, the left kernel is an exponentially decaying function and would be the equivalent of Lapicque's neuron model. On the right, the kernel is an alpha function.\n",
        "\n",
        "In snnTorch, the spike response model is not directly implemented as a filter. Instead, it is recast into a recursive form such that only the previous time step of values are required to calculate the next set of values. This significantly reduces the memory overhead during learning.\n",
        "\n",
        "The filter adopted is unsurprisingly the alpha function on the right animation above, or equivalently a sum of two exponentials. This results in a membrane potential which peaks at some time delay $t_d$ after the input spike. This is often a desirable feature when training networks that rely on spike timing.\n",
        "\n",
        "<center>\n",
        "<img src='https://github.com/jeshraghian/snntorch/blob/master/docs/_static/img/examples/tutorial2/2_9_alpha.png?raw=true' width=\"600\">\n",
        "</center> \n",
        "\n",
        "\n",
        "As the membrane potential is now determined by the sum of two exponentials, each of these exponents has their own independent decay rate. $\\alpha$ defines the decay rate of the positive exponential, and $\\beta$ defines the decay rate of the negative exponential. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3OiADDv7Z_z7"
      },
      "source": [
        "alpha = 0.8\n",
        "beta = 0.7\n",
        "\n",
        "# initialize neuron\n",
        "lif6 = snn.Alpha(alpha=alpha, beta=beta, threshold=0.5)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K8yVKbpyZV4J"
      },
      "source": [
        "Using this neuron is the same as the previous neurons, but the sum of two exponential functions requires the synaptic current `syn` to be split into a `syn_exc` and `syn_inh` component:\n",
        "\n",
        "**Inputs**\n",
        "* `spk_in`: each input voltage spike $S_{\\rm in}[t]$ is sequentially passed in\n",
        "* `syn_exc`: excitatory post-synaptic current $I_{\\rm syn_exc}[t]$ at the present time $t$\n",
        "* `syn_inh`: inhibitory post-synaptic current $I_{\\rm syn_inh}[t]$ at the present time $t$\n",
        "* `mem`: membrane potential $U_{\\rm mem}[t]$ at the present time $t$\n",
        "\n",
        "**Outputs**\n",
        "* `spk_out`: output spike $S_{\\rm out}[t+1]$ at the next time step ('1' if there is a spike; '0' if there is no spike)\n",
        "* `syn_exc`: excitatory post-synaptic $I_{\\rm syn_exc}[t+1]$ at the next time step $t$\n",
        "* `syn_inh`: inhibitory post-synaptic current $I_{\\rm syn_inh}[t+1]$ at the next time step $t$\n",
        "* `mem`: membrane potential $U_{\\rm mem}[t+1]$ at the next time step\n",
        "\n",
        "As with all other neuron models, these must be of type `torch.Tensor`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Bf2fjxBZajeT"
      },
      "source": [
        "# input spike: initial spike, and then period spiking \n",
        "spk_in = (torch.cat((torch.zeros(10), torch.ones(1), torch.zeros(89), (torch.cat((torch.ones(1), torch.zeros(9)),0).repeat(10))), 0) * 0.85).unsqueeze(1)\n",
        "print(f\"spk_in contains {spk_in.size(1)} sample of data across {spk_in.size(0)} time steps.\")\n",
        "\n",
        "# initialize parameters - arg '1' is passed to indicate just one sample of data\n",
        "syn_exc, syn_inh, mem = lif6.init_alpha()\n",
        "mem_rec = []\n",
        "spk_rec = []\n",
        "\n",
        "# run simulation\n",
        "for step in range(num_steps):\n",
        "  spk_out, syn_exc, syn_inh, mem = lif6(spk_in[step], syn_exc, syn_inh, mem)\n",
        "\n",
        "  mem_rec.append(mem.squeeze(0))\n",
        "  spk_rec.append(spk_out.squeeze(0))\n",
        "\n",
        "mem_rec = torch.stack(mem_rec)\n",
        "spk_rec = torch.stack(spk_rec)\n",
        "\n",
        "# Generate Plots\n",
        "fig, ax = plt.subplots(3, figsize=(8,6), sharex=True, \n",
        "                       gridspec_kw = {'height_ratios': [0.4, 1, 0.4]})\n",
        "\n",
        "# Plot input current\n",
        "splt.raster(spk_in, ax[0], s=400, c=\"black\", marker=\"|\")\n",
        "ax[0].set_ylabel(\"Input Spikes\")\n",
        "ax[0].set_title(\"Alpha Neuron Model With Input Spikes\")\n",
        "ax[0].set_yticks([]) \n",
        "\n",
        "# Plot membrane potential\n",
        "ax[1].plot(mem_rec.detach())\n",
        "ax[1].set_ylim([0, 0.6])\n",
        "ax[1].set_ylabel(\"Membrane Potential ($U_{mem}$)\")\n",
        "ax[1].axhline(y=0.5, alpha=0.25, linestyle=\"dashed\", c=\"black\", linewidth=2)\n",
        "plt.xlabel(\"Time step\")\n",
        "\n",
        "# Plot output spike using spikeplot\n",
        "splt.raster(spk_rec, ax[2], s=400, c=\"black\", marker=\"|\")\n",
        "ax[2].set_yticks([])\n",
        "ax[2].set_ylabel(\"Output Spikes\")\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "63V7WsZFa_Uo"
      },
      "source": [
        "As with the Lapicque and Synaptic models, the Alpha model also has options to modify the threshold and reset mechanism."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PYWhh7idiS9t"
      },
      "source": [
        "# Conclusion"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BGCY7g7R0Usy"
      },
      "source": [
        "Now you should understand the basics of several LIF neuron models, how to simulate them, and how to build your own feedforward spiking neural networks.\n",
        "\n",
        "For reference, the documentation [can be found here](https://snntorch.readthedocs.io/en/latest/snntorch.html).\n",
        "\n",
        "In the next tutorial, you will learn how to train these networks to classify spiking and static MNIST datasets."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OTRCyd0Xa-QK"
      },
      "source": [
        "## Further Reading\n",
        "* [snnTorch documentation](https://snntorch.readthedocs.io/en/latest/snntorch.html) of the Lapicque, Leaky, Synaptic, and Alpha models\n",
        "* [*Neuronal Dynamics:\n",
        "From single neurons to networks and models of cognition*](https://neuronaldynamics.epfl.ch/index.html) by\n",
        "Wulfram Gerstner, Werner M. Kistler, Richard Naud and Liam Paninski.\n",
        "* [Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems](https://mitpress.mit.edu/books/theoretical-neuroscience) by Laurence F. Abbott and Peter Dayan"
      ]
    }
  ]
}